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SUMMARY 


The  current  interest  in  aircraft  exhibiting  leading  edge  vortex  flows 
during  takeoff,  landing,  and  high  angle  maneuver  has  led  to  the  computa- 
tional method  development  discussed  in  this  report.  These  methods  are 
directed  to  the  analysis  of  the  forces  and  flow  field  around  arbitrary 
swept  wing  planforms  having  thick  airfoils  with  rounded  leading  edges  and 
leading  edge  vortex  flow.  Due  to  the  nature  of  this  type  of  flow  field 
viscous  effects  must  be  considered  in  the  method.  These  are  not  feasible 
to  obtain  by  solving  the  Navier-Stokes  equations  for  the  entire  flow  field. 
Therefore,  in  contrast  to  other  known  empirical,  potential  flow  methods, 
the  method  discussed  herein  provides  an  optimum  combination  of  viscous  and 
potential  flow  models.  This  combination  is  developed  with  three  goals  in 
mir.d.  One  is  to  maintain  minimum  computer  time,  yet  still  provide  insight 
into  the  mechanisms  of  the  flow  field.  The  second  is  for  the  solutions  to 
be  valid  for  large-scale  configurations,  i.e.  without  Reynolds  number 
limitations.  The  third  goal  is  to  provide  a design  tool  which  by  optimiz- 
ing leading  edge  shape  takes  advantage  of  the  leading  edge  vortex  in 
improving  the  economy,  performance  and  flying  qualities  of  future  aircraft. 

The  method  described  herein  is  a computational  package  which  includes 
four  major  programs  run  with  proper  interfacing  in  an  iterative  cycle. 

These  programs  are:  (1)  a viscous  program  which  solves  parabolized  Navier- 
Stokes  equations  within  a box  containing  the  vortex  flow  on  the  upper  wing 
surface,  (2)  a three-dimensional,  second-order  boundary  layer  program  which 
provides  the  strength  and  distribution  of  the  leading  edge  vorticity  being 
fed  into  the  viscous  box,  (3)  a potential  flow  program  which  models  the 
three-dimensional  thick  wing  using  the  Hess  surface  singularity  method, 
and  (4)  a potential  influence  vortex  program  to  model  the  leading  edge 
vortex  system.  The  two  potential  flow  programs  provide  both  the  pressure 
coefficients  required  for  the  boundary  layer  program  and  the  velocity 
boundary  conditions  for  the  viscous  box.  The  viscous  box  provides  not  only 
the  primary  vortex,  but  other  viscous  effects  such  as  the  vortex  core 
velocities  and  the  secondary  vortex. 
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To  develop  the  interfacing  and  convergence  techniques  of  these  four 
programs,  an  application  was  made  to  a 65°  delta  wing  with  spanwise  varia- 
tion in  leading  edge  radius.  Operating  experience  is  presented  in  this 
report.  Each  cycle  through  all  the  programs  has  contributed  not  only  to 
the  operating  techniques  for  interfacing  the  programs,  but  to  a knowledge 
of  the  critical  parameters.  It  is  shown  in  this  report  that  aircraft  loads 
and  moments  can  be  critically  dependent  on  the  strength  and  distribution  of 
the  vorticity  shed  from  the  wing  leading  edge  into  the  primary  vortex. 

This  vorticity  is  dependent  on  local  angle  of  attack,  leading  edge  shape, 
Reynolds  number,  and  the  boundary  layer  state  (laminar  or  turbulent). 


Several  trial  runs  have  been  made  exercising  all  program  elements 
sequentially  to  develop  effective  interface  and  convergence  techniques. 
Finally,  these  techniques  are  employed  in  a case  where  three  consistent 
cycles  of  the  full  program  are  completed  to  evaluate  the  extent  of  solution 


convergence.  This  case  suggests  additional  convergence  techniques  while 
demonstrating  the  feasibility  of  the  overall  method  and  providing  direction 
for  future  work. 


Recommendations  presented  for  future  efforts  lie  in  two  categories. 

One  concerns  program  changes  to  make  the  method  more  versatile  and  accurate. 
The  other  reflects  on  the  need  for  both  additional  operational  experience 
with  the  current  program  and  experimental  verification  of  some  of  the 
resulting  flow  field  parameters. 
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Lift  coefficient 
P i tch i ng  moment 
Pressure  coefficient 
Wing  root  chord 

Momentum  coefficient  for  spanwise  blowing 
Mach  number 

Direction  normal  to  surface 
Pressure 

Radius  of  free  vortex 
Reynolds  number 

Distance  along  stream  tube  in  boundary  layer 

velocities  in  x,  y,  z direction,  respectively 

local  edge  velocity  of  boundary  layer 

component  of  free  stream  normal  to  wing  leading  edge 

distance  normal  to  wing  leading  edge 

rectangular  coordinates  - depends  on  program. 

See  Figure  2.1  for  Hess  program  dnd  viscous  box. 

See  Figure  4.2  for  boundary  layer. 

intersection  of  box  plane  with  root  chord 

c i rcu I at  ion 

boundary  layer  thickness 

angle  between  box  no-slip  surface  and  airfoil  chord 
wing  angle  of  attack 

density  and  radius  of  airfoil  leading  edge 
kinematic  viscosity 

parameter  defining  differencing  technique 


1 1 


Subscripts 


wing  sweep 

fraction  of  spanwi se  position  and  vorticity  normal 
to  viscous  box  x-z  plane 

slope  of  ai rfoi 1 

vorticity  normal  to  y-z  plane  of  viscous  box 
vorticity  normal  to  x-y  plane  of  viscous  box 
total  vorticity 


related  to  free  vortex 
wing  thickness 


1.0  INTRODUCTION 


The  favorable  lifting  effects  of  leading  edge  vortices  on  delta  wings 
have  been  under  investigation  for  many  years.  Several  years  ago  the  TFX 
configurations  resulting  in  the  Fill  found  favorable  effects  of  the  leading 
edge  "strake"  effect  on  maximum  lift  coefficients.  The  FI  1 6 , F14  and  other 
fighter  type  of  aircraft  have  more  recently  taken  advantage  of  vortex  lift. 
Latest  application  of  vortex  lift  is  found  in  the  "Space  Shuttle"  configu- 
ration and  the  supersonic  "Concorde".  Future  supersonic  configurations 
will  be  aided  or  plagued  by  the  effects  of  leading  edge  vortex  formation 
depending  on  the  point  of  view  or  application  thereof.  Vortex  flow 
contains  considerable  energy  which  most  often  shows  up  as  drag;  however, 
the  current  direction  is  to  use  vortex  flow  to  improve  lift.  The  question 
remaining  is,  "Can  L/D  also  be  improved  when  designing  for  vortex  lift?". 
For  many  cases  the  leading  edge  vortex  is  persistent  regardless  of  attempts 
to  eliminate  it;  therefore,  it  must  be  put  to  useful  work. 

There  have  been  many  attempts  to  model  the  flow  for  delta  wings  with 

leading  edge  vortices.  Matoi^  has  provided  a good  survey  of  many  of  these 

attempts.  The  efforts  of  most  have  been  limited  to  potential  flow  models, 

assuming  rotational  elements  for  viscous  effects  such  as  the  leading  edge 

feeding  sheet.  The  basic  assumption  of  slender-body,  conical  flow  have 

been  followed  by  most  investigators  such  as  Mangier  and  Smith^.  Results 

of  these  early  investigators  have  laid  a good  foundation  for  further 

research,  but  they  are  not  successful  in  obtaining  forces,  and  moments 

o U 

over  a general  range  of  configurations.  Polhamus^  and  Lamar  have  had 
good  success  in  predicting  vortex  lift  with  the  leading  edge  suction 
analogy  for  sharp  edge  wings.  This  method,  however,  does  not  describe  the 
pressures  or  flow  mechanism.  Good  results  were  obtained  most  recently  by 
Weber,  et  al.'’  and  Rao  and  Nathnan^,  where  the  trailing  edge  conditions  of 
delta  wings  are  included  in  the  boundary  conditions.  Kandil?  has  also 
provided  a good  potential  flow  model  and  convergence  criteria.  Currently 
these  programs  appear  to  have  limits  that  restrain  their  application  to 
certain  configurations.  These  limits  are  as  follows: 
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(1)  Thin  wings  only  are  considered. 

(2)  Starting  conditions  at  leading  edge  must  be  assumed. 

(3)  The  viscous  condition  of  no  slip  at  the  wing  surface 
in  the  area  of  the  free  vortex  is  not  satisfied  and 
this  can  affect  final  vortex  position  and  strength. 

(4)  No  secondary  vortex  system  is  included. 

(5)  Vortex  burst  cannot  be  predicted. 

Based  on  verbal  comments  from  others  involved  in  this  field,  some  of  the 
above  investigators  are  beginning  to  apply  the  necessary  viscous  tools  to 
remove  the  above  limitations.  G.  F.  Hall,  et  al.®  have  recognized  the 
need  to  combine  viscous  solutions  with  potential  flow  techniques  for 
application  to  the  formation  of  a tip  vortex  from  a semi-infinite  wing. 
Their  methods  are  similar  to  those  employed  by  Lockheed.  Lockheed-Georg ia 
Company  has  been  developing  the  necessary  viscous  tools  for  approximately 
three  years  under  Office  of  Nava)  Research  contract  and  now  has  integrated 
the  computer  programs  that  have  the  potential  to  remove  practically  all  of 
the  above  limitations.  The  fifth  limitation,  vortex  bursting,  is  not  well 
resolved,  but  the  methods  do  give  indication  of  incipient  burst.  Some 
empirical  work  must  still  be  applied  to  this  prediction.  The  resulting 
hybrid  viscous/potential  flow  program  is  intended  to  provide  the  flow 
field  velocities,  pressures  and  forces  for  arbitrary  configurations  having 
interfering  vortex/wing  flows.  Thick  rounded  leading  edge  wings  with 
leading  edge  vortex  due  to  flow  separation  from  the  swept  leading  edges 
can  be  analyzed.  This  hybrid  program  combines  two  elliptic  potential 
flow  programs,  a "viscous  box"  Navi er-Stokes  solution  containing  the 
vortex  and  a boundary  layer  program. 

The  original  effort  to  develop  the  "viscous  box"  technique  is 
reported  in  references  9 and  10.  These  reports  were  primarily  concerned 
with  leading  edge  vortex  control  by  spanwi se  blowing^  SWB.  Application 
was  limited  to  straight  or  moderately  swept  wings  since  the  box  boundary 
conditions  and  flow-field  velocities  were  obtained  by  a Green's  function 
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integral  technique.  Later  the  interest  turned  to  highly  swept  wings 
without  SWB.  This  required  a more  accurate  technique  to  supply  the 
boundary  conditions  to  the  viscous  box.  Reference  11  reports  the  develop- 
ment of  the  hybrid  viscous  box/elliptical  potential  flow  vortex  lattice 
method  as  applied  to  thin  sharp  edged  delta  wings.  In  this  technique  the 
elliptic  vortex  lattice  program,  through  iteration,  supplied  the  box 
boundary  velocities,  but  the  vorticity  shed  from  the  leading  edge  across 
the  box  boundary  still  had  to  be  empirically  estimated.  A sensitivity 
study  showed  the  strength  and  direction  of  this  feeding  vorticity  to  be 
very  critical  to  the  overall  solution.  This  led  to  the  development  of  a 
boundary  layer  program  to  supply  this  need.  With  this  additional  tool, 
thick  rounded  leading  edge  wings  with  leading  edge  separation  can  be 
treated  as  well  as  thin  wings  if  joined  with  a thick  wing  potential  flow 
model  and  the  viscous  box  concept.  The  development  and  integration  of 
these  tools  for  the  thick  wing  is  the  subject  of  this  report. 


2.  SOLUTION  PROCEDURE 


The  computational  package  includes  four  major  programs  which  are  run 
with  proper  interfacing  in  an  iterative  cycle.  The  coordinate  system  for 
these  programs  and  the  regions  where  they  are  used  are  indicated  in  figure 
2.1  and  2.2.  The  flow  chart  of  figure  2.3  shows  the  sequence  of  computa- 
tions using  these  programs.  The  capabilities  of  individual  programs  are 
described  in  subsequent  chapters;  the  overall  computational  sequence  and 
the  interfacing  computational  work  between  the  major  calculations  are 
described  below. 

The  major  programs  employed  in  the  computations  are:  (l)  Aircraft 
Interference  Program  (AIP),  which  models  the  three-dimensional  thick  wing 
using  potential  flow  assumption,  (2)  a vortex  lattice  program  to  model  the 
leading  edge  vortex  system  for  the  delta  wing,  (3)  a three-dimensional, 
second-order,  boundary  layer  program,  and  (A)  a viscous  box  program  that 
solves  parabolized  Navier-Stokes  equations  on  the  upper  surface  of  the 
delta  wing.  The  potential  flow  models  provide  the  pressure  coefficients 
required  for  the  boundary  layer  calculation  in  addition  to  the  velocity 
boundary  conditions  for  the  viscous  box  calculations.  The  step-by-step 
procedure  given  below  supplements  the  flow  chart  of  figure  2.3. 

(1)  Define  coordinates  of  panels  on  the  wing  for  AIP  program. 
Determine  the  centroids  of  these  panels  which  form  the  control  points  for 
AIP  program. 

(2)  Define  control  points  for  the  vortex  lattice  program.  An 
optimum  choice  of  panelling  procedure  for  these  two  steps  is  described  in 
Section  3. 

(3)  Define  in  global  coordinates,  for  AIP  program,  the  coordinates 
of  grid  points  on  the  box  boundaries. 
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Figure  2.2.  Airfoil/Viscous  Box  Relation 


(A)  Define  leading  edge  vortex  strength  and  position  for  starting 
the  computations  at  zeroth  iteration.  This  is  estimated  based  on  experi- 
mental or  other  available  data  (see  Section  3).  For  subsequent  iteration 
this  is  defined  using  the  final  results  of  Nav ier-Stokes  computations. 

(5)  Run  Vortex  Lattice  Program  and  obtain  surface  velocities  at 
lattice  control  points  defined  in  Step  2. 

(6)  Run  an  interfacing  program  and  using  the  data  of  Step  5 above, 
obtain  velocities  due  to  leading  edge  vortex  system  at  AIP  control  points 
def i ned  i n Step  1 . 

(7)  Run  AIP  using  data  of  Step  1 with  the  velocities  obtained  in 
Step  6 imposed  as  additive  boundary  condition  on  the  surface.  Obtain 
surface  pressure  coefficients. 

(8)  To  determine  the  velocity  boundary  condition  for  the  viscous 
box  solution,  run  Vortex  Lattice  Program  first  to  compute  velocities  at 
the  box  grid  points  defined  in  Step  3 and  save  the  data.  Then  run  the 
AIP  program  which  is  designed  to  add  the  effect  of  leading  edge  vortex  to 
the  velocities  induced  by  the  surface  singularities. 

(9)  Run  an  interfacing  program  to  generate  surface  ordinates  and 
pressure  distribution  along  lines  normal  to  the  wing  leading  edge  for  the 
box 'solution  planes. 

(10)  Using  data  of  Step  9 above,  run  the  boundary  layer  program  for 
the  viscous  solution  planes.  Obtain  the  leading  edge  vorticity  that  is 
fed  onto  the  upper  surface  of  the  wing  from  these  runs  at  each  box 
solution  plane. 

(11)  The  first  box,  run  the  viscous  box  program  at  initial  plane. 

Use  very  small  AX  values  and  large  number  of  iterations,  approximately  20, 
between  velocity  transport  equation  and  Poisson's  equation  for  velocity. 
Define  the  resultant  solution  as  the  initial  conditions  (see  Section  5-2). 
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(12)  With  the  initial  conditions  obtained  above  in  Step  11,  the 
velocity  boundary  conditions  determined  in  Step  8,  and  the  leading  edge 
vorticity  values  of  Step  10,  run  the  viscous  box  program  for  all  the 
boxes. 

(13)  From  the  vorticity  distribution  of  viscous  solution,  define  a 
new  vortex  position  and  strength.  If  the  position  and  strength  have  not 
converged  within  acceptable  limits,  repeat  calculations  from  Step  k.  If 
they  have  converged,  run  the  potential  models  alone  once  to  obtain  final 
pressure  distribution  on  the  surface  and  the  force  and  moment 
coefficients. 
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3.  POTENTIAL  FLOW  MODELS 

The  potential  flow  calculations  are  carried  out  using  two  separate 
programs.  A vortex  lattice  program  models  the  leading  edge  vortex  system 
and  computes  the  velocities  induced  by  such  a system  at  chosen  points  in 
the  field.  The  velocities  computed  on  the  wing  surface  are  imposed  as 
additional  boundary  condition  in  a second  program  known  as  Aircraft 
Interference  Program  (AIP),  which  models  the  thick  wing  using  singularity 
distributions  on  the  wing  surface. 


Figure  3-1  illustrates  the  free  vortex  system.  For  clarity  only  one 
vortex  element  is  shown  while  the  actual  computations  use  20  such  vortex 
elements.  The  bound  part  of  the  vortex  filament  is  located  at  the  wing 
center  plane  and  the  free  element  is  allowed  to  change  direction  when  it 
meets  the  other  filaments.  To  start  with  the  vortex  geometry  and  the 
circulation  strengths  are  estimated  from  data  available  through  experiments 
or  any  other  empirical  procedure.  On  subsequent  cycles  of  iteration,  the 
vortex  system  is  reconstructed  from  the  vorticity  distributions  obtained  by 
the  viscous  box  solution.  Once  the  vortex  geometry  and  the  circulation 
strengths  are  completely  specified,  the  vortex  lattice  program  computes 
the  velocity  at  any  chosen  point  in  the  field  using  Bio-Savart's  law. 

The  initial  free  vortex  core  position  and  strength  has  been  estimated 
by  the  methods  described  in  reference  11  with  an  additional  modification 
to  account  for  the  effect  of  a rounded  leading  edge.  The  effect  of  the 
rounded  leading  edge  is  obtained  from  evaluation  of  experimental  data 
comparing  sharp  and  rounded  leading  edge  in  reference  12.  Resulting 
equations  for  these  vortex  characteristics  are: 


0.5380  i . 0.01281  <£>*  - 0.007226  (i|’  - 0.1268  l±) 

where  Yy  is  the  horizontal  position 

Z\j  is  the  vertical  position 

Zt  is  the  local  wing  semi-thickness 

Cr  is  the  root  chord 

X is  distance  along  root  chord  from  apex. 

The  vortex  core  radius  has  been  estimated  based  on  the  laser  data  obtained 
from  previous  work  discussed  in  reference  11.  Assuming  the  core  radius  is 
proportional  to  T of  the  vortex  or  sin2otcosa,  the  core  size  is  estimated 
for  a = l!*°  as  follows. 

R V . x 

= 0.02008  -jr- 
CR  CR 

This  value  helps  to  determine  the  viscous  box  size. 

3.2  Thick  Wing  Potential  Flow  Model 

The  potential  flow  past  the  thick  wing  is  analyzed  using  a program 
referred  to  as  the  Aircraft  Interference  Program  (AIP)  which  is  a develop- 
ment of  the  method  originally  formulated  by  Hess  (reference  13)  under  Navy 
sponsorship.  This  program  uses  combinations  of  source  and  doublet  panels 
to  compute  subsonic  inviscid  flows  about  arbitrary  configurations.  The 
Kutta  condition  on  lifting  components  is  satisfied  by  requiring  the  upper 
and  lower  surface  pressures  to  be  equal  near  the  trailing  edge.  This  has 
been  found  to  yield  better  results  than  the  specification  of  trailing  wake 
posi tion. 


Among  the  major  modifications  to  the  original  Hess  program  that  have 
been  made  at  Lockheed-Georg ia  are: 


(1)  Improved  coding  to  eliminate  errors  in  the  originally 
acquired  program  and  to  make  the  program  computationally 
more  efficient. 

(2)  Modified  program  input  and  output  making  it  easier  to  use. 

(3)  Development  of  a configuration  loft  program  to  automatically 
determine  component  intersections  and  to  permit  visual 
verification  of  the  paneling. 

The  AIP  program  has  the  capability  of  providing  surface  velocity 
components  and  pressures,  aerodynamic  forces  and  moments,  and  off-body 
flow  field  information.  The  program  is  currently  applicable  to  arbitrary 
configurations  for  which  viscous  effects  can  be  neglected.  Figures  3.2 
and  3*3  show  a typical  transport  configuration  and  attendant  wing  surface 
pressures.  The  lack  of  agreement  between  theory  and  experiment  evident 
in  the  lower  surface  pressure  comparisons  occurs  because  the  pylons  and 
nacelles  were  not  included  in  the  theoretical  calculations.  The  slight 
over-prediction  of  upper  surface  pressures  is  attributable  to  the 
neglect  of  viscous  effects.  The  good  agreement  between  theory  and 
experiment  evident  in  the  vicinity  of  the  leading  edge  is  of  particular 
significance  for  the  subject  study.  These  calculations  were  performed 
with  528  surface  panels  and  require  approximately  1 6 1 seconds  computation 
time  on  a CDC  7600  computer.  Less  computer  time  is  required  for  lower 
Mach  numbers. 

3.3  Integration  of  Vortex  Lattice  Model  and  AIP  Program 

The  effect  of  the  leading  edge  vortex  system  on  the  wing  is 
established  by  the  inclusion  of  its  induced  velocities,  as  determined  by 
the  Lattice  Program,  in  the  boundary  conditions  at  the  AIP  program's 
collocation  points  on  the  wing.  The  computed  results  of  AIP  program  will 
then  reflect  the  presence  of  the  leading  edge  vortex.  Also,  at  the 
boundary  grid  points  of  viscous  box  program,  the  velocities  are  determined 
by  adding  the  velocities  at  these  points  obtained  through  AIP  program  as 
stated  above  to  the  velocities  induced  by  the  leading  edge  vortex  system 
i tsel f . 
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THEORY 


Figue  3-3  Pressure  Distributions  for  C-5A  Wing  at  36.9%  Span 


WfW— ■ ' MSI 

j9i 

B<?£«PtE3 

s 

V*|R 

Initially  the  two  potential  flow  programs  were  organized  to  run 
separately  with  an  additional  peripheral  program  to  integrate  the  results. 
This  was  considered  a desirable  feature  during  the  development  stage,  but 
these  two  codings  were  combined  into  one  for  case  4 runs  discussed  in 
Section  6.3.  A diagrammatic  representation  of  the  coupling  process 
currently  in  use  is  given  in  figure  3-^- 

One  major  aspect  of  the  interfacing  procedure  between  these  two 
programs  needs  to  be  mentioned.  The  conventional  method  of  paneling  wings 
and  other  lifting  bodies  in  the  AIR  program  is  by  means  of  "n  lines"  which 
are  lines  drawn  on  the  surface  and  oriented  roughly  in  the  direction  of 
the  flow.  The  space  between  two  of  these  n lines  is  called  a strip  which 
is  divided  in  the  chordwise  direction  to  form  panels.  All  calculations 
are  then  carried  out  at  control  points  which  are  situated  at  the  centroids 
of  these  panels.  A conventional  paneling  laid  out  on  a delta  wing  is  as 
shown  in  figure  3-5.  This  configuration  is  not  entirely  compatible  with 
the  vortex  program  since  the  bound  vortices  would  pass  very  close  to  the 
control  points  giving  uneven  velocity  distribution.  On  the  other  hand, 
a paneling  with  the  conical  flow  n-lines  converging  to  the  apex  of  the 
wing,  as  shown  in  figure  3.6,  for  the  A I P program,  was  found  to  give 
uneven  variations  of  the  field  variables.  Further,  it  also  proved  very 
sensitive  to  the  paneling;  i.e.  a small  change  in  the  panel  sizes  at  the 
apex  caused  a 100%  change  in  lift  coefficient. 

So,  it  was  decided  to  revert  to  conventional  panelling  and  after  some 
experimentation,  a configuration  with  218  panels  was  chosen  and  provided 
satisfactory  accuracy  without  undue  computational  time.  To  prevent 
interference  between  vortex  elements  of  the  thick  wing  model  and  the  free 
vortex  model,  the  free  vortex  model  is  run  independently  with  its  own 
optimum  panels  and  velocities  computed  at  their  centroids.  The  surface 
velocities  at  the  centroids  of  the  thick  wing  panels  due  to  the  free 
vortex,  are  obtained  by  interpolation  and  used  in  the  boundary  condition 
for  the  thick  wing  model. 
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A.I.P.  EXECUTED 
PARTIALLY 


OUTPUT: 

ON-BODY  LOCATIONS 
STORED  IN  ONBODPTS 


VORTEX  PROG  RUN 
FOR  ON-BODY  PTS 


OUTPUT: 

ON-BODY  VELOCITIES 
STORED  AS  VORONBOOVEL 


OFF-BODY  POINTS 
GEOMETRY  PROGRAM 

OUTPUT: 

OFF-BODY 
STORED  AS 

LOCATIONS 

OFFBODPTS 

VORTEX  PROG  RUN 
FOR  OFF-BODY  PTS 


OUTPUT: 

OFF-BODY  VELOCITIES 
STORED  AS  VOROFFBODVEL 


VELOCITY  COMBINATION 
PROGRAM 


AIP 

INPUT: 


(a)  WING  GEOMETRY 

(b)  VELOCITIES  DUE 
TO  VORTEX  ON-  5 
OFF-BOOY 

(c)  OFF-BODY  PTS 


OUTPUT 


OFF-BODY  VELOCITIES  -«■ 
FOR  VORTEX  BOX  PROGRAM 


ON-BODY  CPS  TO  BE  USED  BY 
BOUNDARY  LAYER  PROGRAM 


Figure  3.^  Potential  Flow  Integration  Program 
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4.0  BOUNDARY  LAYER  PROGRAM 


4.1  Requirement  for  Leading  Edge  Vorticity 

As  indicated  in  reference  11  and  in  Sections  1.0  and  2.0,  it  is 
imperative  to  know  the  strength  and  spatial  orientation  of  the  vorticity 
shed  at  the  leading  edge.  This  parameter  is  a function  of  local  angle  of 
attack,  circulation  and  the  shape  of  the  leading  edge.  These  conditions 
in  turn  determine  the  local  pressure  distribution  and  the  boundary  layer 
characteristics.  Since  the  configurations  of  concern  are  usually  swept 
and  three-dimensional,  a 3"0  boundary  layer  program  that  can  predict  the 
boundary  layer  characteristics  very  near  flow  separation  is  required.  In 
some  cases  vortex  shedding  (flow  separation)  does  not  occur  until  slightly 
aft  of  the  leading  edge  (for  airfoils  having  leading  edge  separation  as 
opposed  to  trailing  edge  separation).  For  these  cases  it  is  useful  and 
more  accurate  to  have  the  boundary  layer  as  thick  as  possible  in  order  to 
cross  a multiple  number  of  viscous  box  mesh  points  as  the  boundary  layer 
enters  the  box.  Strictly  speaking,  it  is  required  to  compute  a boundary 
"region",  i.e.  a region  including  higher-order  effects  than  those 
occurring  in  first-order  boundary  layer  theory  (see  figure  4.1).  In  order 
to  establish  with  accuracy  the  forward  edge  velocity  and  vorticity  values 
for  the  box,  it  is  necessary  to  perform  a computation  from  the  wing 
attachment  line  through  the  boundary  layer  and  into  the  second  order, 
rapidly  thickening,  boundary  layer  region.  This  technique  not  only 
provides  a thick  boundary  layer  but  it  brings  the  boundary  layer  solution 
very  close  to  separation.  It  is  assumed,  therefore,  when  the  forward 
marching,  parabolic  boundary  layer  solution  "blows-up"  the  boundary  layer 
separates.  This  allows  one  to  determine  at  high  Reynolds  number  if  there 
is  in  fact  leading  edge  separation.  Recognizing  the  sensitivity  of  the 
viscous  box  vortex  flow  region  to  the  vorticity  imposed  by  the  leading 
edge  boundary  layer  region,  it  is  essential  that  three  dimensionality  be 
accounted  for  in  any  model  of  the  boundary  region.  For  this  purpose,  it 
is  deemed  necessary  and  sufficient  to  apply  the  "infinite  yawed  wing" 
assumption.  This  method  has  proven  quite  effective  in  the  analysis  of 
boundary  layers  near  attachment  lines  and  even  for  some  complete  wings  of 
sufficiently  high  aspect  ratio,  as  shown  by  Nash  and  Tseng1^.  Using  the 
three-dimensional  boundary  layer  methods  of  Nash  and  Patel '5,  Nash  and 
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2nd  ORDER 
REGION 


Figure  4.1  Section  Showing  Boundary  Layer  First  and 
Second  Order  Domains 


Figure  4.2  Wing  Section  Showing  Boundary  Region  Coordinate  System 


Scruggs  as  consultants  to  Lockheed  developed  the  infinite  yawed  wing, 
second-order  boundary  layeif  program  for  the  current  effort. 


4.2  Second-Order  Boundary  Layer  Method 

In  order  to  approximate  the  flow  upstream  of  the  viscous  box 
boundary,  the  usual  boundary  layer  assumptions  must  be  augmented  by  the 
appropriate  second-order  effects.  Since  the  boundary  region  experiences 
rapid  increase  in  thickness  and  considerable  turning  of  the  flow  in  the 
spanwise  direction,  two  additional  effects  must  be  represented:  three- 
dimensionality  and  normal  pressure  gradient.  For  the  first  of  these 
requirements  it  was  deemed  adequate  to  assume,  at  least  locally,  that  the 
flow  is  of  the  infinite  yawed  cylinder  type.  This  retains  three- 
dimensionality  but  permits  much  faster  computation  of  the  flow  field.  In 
this  case  the  metric  on  the  wing  surface  is 

ds2  = (hjdx)2  + (h2dy)2  + dz2  (4.1) 

where,  as  shown  in  figure  4.2  the  coordinate  x is  measured  along  lines 
normal  to  the  wing  leading  edge  and  z is  measured  along  the  generators. 

The  first-order  boundary  layer  momentum  equations  are  then 


JLi£  + JLiy.+  _L_L(£)  + _L  _ _v_  a2u 
hj  3x  h2  3y  hj  3x  p h2  3y  ^2  jy2 


(*.2) 


U 3w  _V_  _3w  _1_  3wv  _ _v_  32w 

hi  3x  h2  3y  h2  3y  h \ 9Y 


The  turbulent  kinetic  energy  equation  is  also  solved  to  obtain  the 
turbulent  shear  stresses.  The  we  1 1 -documented  method  of  Nash1-*  is  used 
to  model  the  shear  stresses  in  terms  of  turbulent  kinetic  energy.  The 
continuity  equation  is  used  to  solve  for  normal  velocity,  V. 

In  order  to  allow  for  normal  pressure  gradients  in  the  case  of  a 
rapidly  thickening  domain,  a second  order  approximation  is  required. 

This  takes  the  form  of  an  approximation  to  the  normal  momentum  equation: 
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k21uv  - k12u*  + ^(f)  = 0,  (4.4) 

where  K2j  and  Kj2  are  curvatures  of  the  coordinate  system.  The  numerical 
method  utilizes  a curvilinear  coordinate  system,  which  is  evolved  as  the 
flow  develops  in  the  downstream  direction.  Thus  (see  figure  4.1)  the  two 
curvatures  can  be  computed  at  each  step  of  the  solution  in  the  stream 
direction.  The  continuity  equation  in  this  system  takes  the  form, 


J_  3U  + J_  3V 

hj  3x  h2  3y 


k2!u 


0. 


(4.5) 


The  procedure  used  to  solve  the  second  order  approximation  involves  a 
retrospective  treatment  of  equation  (4.4).  The  basic  scheme  used  for  the 
first-order  boundary  layer  solution  is  well  tested  and  in  this  manner  the 
essential  features  of  that  scheme  can  be  retained.  The  schematic  in 
figure  4.3  summarizes  the  iterative  procedure  used  to  compute  the  higher 
order  flow.  Pressures  provided  by  the  potential  flow  solution  are  applied 
at  the  wing  surface.  Following  a first  solution  of  the  boundary  layer 
equations,  an  incremental  pressure  gradient  is  found  retrospectively  from 
the  normal  momentum  equation  such  that  the  flow  curvatures  are  balanced. 
The  incremental  pressure  variation  through  the  layer  is  found  by  inte- 
grating the  normal  momentum  equation  from  the  wall  to  the  outer  edge  of 
the  layer.  This  increment  is  then  differentiated  with  respect  to  x and 
added  to  the  previous  pressure  gradient  used  in  the  boundary  layer 
equations.  Thus,  3p/3x  appearing  in  equation  (4.2)  becomes  a function  of 
both  x and  y.  The  new  3p/3x  implied  at  the  edge  of  the  layer  is  used 
in  the  Euler  equations  to  obtain  new  edge  values  of  velocity.  The  calcu- 
lation then  returns  to  the  boundary  layer  equations  for  a new  iteration 
and  the  procedure  is  continued  to  convergence. 


The  basic  numerical  procedure  is  similar  to  that  used  in  the  vorticity 
box  method  (Section  5)  but  with  some  modifications  to  treat  the  present 
problem.  The  scheme  is  forward  marching  and  as  such  requires  initial  con- 
ditions. These  initial  conditions  at  the  attachment  line  depend  on  the 
value  of  wing  sweep.  Beyond  about  25°  the  flow  is  turbulent.  For  the  case 
of  a laminar  leading  edge,  exact  attachment  line  solutions  are  available 
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to  start  the  computation.  These  solutions  for  U and  W as  a function  of 
y,  and  near  the  attachment  line,  are  shown  in  figure  4.4.  For  a turbulent 
leading  edge  appropriate  modifications  should  be  made  to  the  profiles. 


It  has  been  shown  experimentally  that  flow  along  the  attachment  line 
of  an  infinite  swept  wing  will  be  turbulent  when  a parameter 

C*  = w~2/v3U_/9x,  1 

e e 

exceeds  a value  of  about  1.4 x 10s  (reference  16).  We  is  the  edge  velocity 
along  the  attachment  line  and  U_  is  the  edge  velocity  normal  to  the  attach- 
ment  line.  Turbulent  flow  was  not  indicated  for  cases  1 through  3 of 
these  investigations  due  to  an  error  in  the  program,  but  case  4 does  predict 
turbulent  flow  as  expected  for  a 65°  wing.  When  turbulence  is  indicated,  it 
is  current  practice  to  approximate  the  attachment  line  profile  via  a power 
law.  For  the  present  application  this  approximation  may  be  too  crude.  If 
computational  experiments  so  indicate,  the  experimental  profiles  of  Cumpsty 
and  Head  (reference  16)  should  be  curve  fitted  and  used  instead. 

4.3  Boundary  Layer  Results 

Calculations  by  this  method  have,  thus  far,  been  well-behaved;  the 
number  of  iterations  required  has  not  increased  and  no  stability  problems 
have  been  encountered.  Figure  4.5  shows  some  typical  results  from  the 
boundary  layer  program  for  the  airfoil  at  the  face  of  box  2 on  the  65° 
delta  wing.  The  calculated  boundary  layer  thickness  is  shown  as  a function 
of  the  airfoil  chord  and  as  plotted  on  the  airfoil  nose.  The  pressure 
gradient  is  highly  adverse  from  the  airfoil  nose  and  aft  as  the  plot  of  Cp 
shows.  Also  illustrated  are  the  points  where  the  flow  is  caused  to  be 
turbulent  and  where  the  normal  momentum  equation  is  put  into  action 
(dCp/dN  0).  Between  X/C  = .0012  and  .0016  the  flow  separates  since  the 
boundary  layer  program  "breaks  down"  at  X/C  = .0016.  A finer  X/C  choice 
aft  of  X/C  “.0012  would  have  been  desirable  to  obtain  the  separation  point 
more  accurately,  and  therefore  be  a little  more  accurate  on  the  rate  of 
vorticity  shed  from  the  boundary  layer.  Figure  4.6  shows  the  boundary 
layer  profile  at  X/C  “.0012,  just  prior  to  separation.  The  inflection 
in  the  profile  indicates  an  approach  to  separation. 
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Figure  4.5.  Upper  Surface  Boundary  Layer  for  Box  2 Face 


More  detailed  investigations  of  the  boundary  layer  program  results 
were  made  on  another  wing  for  which  experimental  pressure  data  are 
available  and  a leading  edge  vortex  is  formed.  This  wing  has  a very  blunt 
leading  edge  but  there  is  leading  edge  separation  making  it  ideal  to  test 
the  boundary  layer  program.  Figure  4.7  shows  the  airfoil  and  experimental 
pressure  distributions. 

The  pressure  distribution  with  the  vortex  effect  is  all  experimental; 
however,  the  adverse  pressure  distribution  for  fully  attached  flow  is 
estimated  by  extrapolation  from  attached  flow  results  at  lower  angles  of 
attack.  The  airfoil  is  the  streamwise  section  of  the  wing  at  36.4%  of  the 
wing  span.  The  wing  is  of  aspect  ratio  4.0,  22°  leading  edge  sweep  and 
0.7  taper  ratio.  The  airfoil  type  is  constant  with  no  twist.  A steady 
leading  edge  vortex  is  formed  through  the  assistance  of  spanwise  blowing 
from  the  wing  root  at  33%  chord.  Details  of  this  test  and  configuration 
are  given  in  reference  17- 

Boundary  layer  analyses  of  this  wing  section  is  accomplished  by  first 
determining  the  section  coordinates  normal  to  the  leading  edge  at  the  36.4% 
span  position  then  assuming  that  the  pressure  isobars  are  parallel  to  the 
constant  X/C  values.  The  infinite  yawed  wing  boundary  layer  analysis  is 
then  applied.  Figure  4.8  shows  the  boundary  layer  thickness  as  a function 
of  the  X^/C  normal  to  the  leading  edge  (C  is  still  the  streamwise  chord; 
therefore,  X^/C  = (X/C)cos  22°).  Two  Reynolds  number  results  are  shown, 
and  for  the  small  change  no  significant  effect  of  Reynolds  number  is  found 
except  the  relative  boundary  layer  thickness.  All  curves  of  figure  4.8 
are  for  laminar  flow,  which  is  reasonable  at  these  Reynolds  numbers.  Two 
of  the  curves  have  the  normal  momentum  equation  operative  after  station  9 
(Xn/C«.023). 


The  low  Reynolds  number  of  figure  4.8  is  the  actual  experimental 
value  based  on  the  streamwise  chord  length  of  8.2834  cm  (0.44  ft).  For 
this  case,  boundary  layer  stations  have  been  analyzed  to  come  close  to 
the  separation  point.  Initially,  no  data  were  obtained  between  stations  9 
and  II  and  this  produced  a slightly  different  curve  (not  shown).  Also, 
well  beyond  station  13  convergence  could  not  be  met;  therefore,  closer 
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Figure  ^ . 8 Boundary  Layer  Thickness  for  22°  Wing 


1 


data  points  around  station  13  were  input  to  the  program.  In  the  second 
trial,  with  station  10,  included  along  with  points  beyond  13,  separation 
is  not  indicated  until  just  slightly  past  X^/C  * .05**.  It  is  interesting  to 
note  that  in  this  case  figure  4.8  shows  that  without  the  normal  momentum 
equation  operating  separation  occurs  at  about  the  same  place.  The  boundary 
layer  thickness  is  somewhat  smaller  though. 

Figure  4.9  shows  the  boundary  layer  profiles  for  stations  9.  14,  15, 
and  16  for  the  lower  Reynolds  number  data.  The  progression  towards  a 
separated  profile  is  obvious  with  that  at  station  16  coming  very  close  to 
the  separation  point.  Ideally  the  program  should  be  run  at  more  points 
until  the  initial  separation  profile  is  found,  i.e.  with  3u/3y=0  at  the 
surface. 

The  separation  point  is  approximated  at  XN/C  = .054  or  X/C  = .058.  As 
shown  on  figure  4.7,  this  is  in  good  agreement  with  the  separation  indi- 
cated by  the  experimental  pressure  distribution.  It  is  important  to  note 
that,  if  the  actual  experimental  pressure  distribution  is  used  with  its 
relatively  low  adverse  pressure  gradient,  no  separation  occurs  in  the 
boundary  layer  program.  It  was  therefore  necessary  to  estimate  the  dis- 
tribution as  shown  in  figure  4.7  for  attached  flow.  After  some  thought, 
this  is  considered  the  most  accurate  technique.  The  boundary  layer 
program  is  designed  for  attached  flow  only,  even  though  the  second-order 
modification  brings  it  close  to  separation;  therefore,  the  correct  tech- 
nique is  exactly  what  is  done  in  the  theoretical  approaches  of  the 

i 

following  sections.  The  attached  flow  pressure  distribution  at  the 
appropriate  angle  of  attack  and  in  the  presence  of  the  leading  edge  vortex 
is  input  to  the  boundary  layer  program.  No  leading  edge  vortex  sheet  is 
included  for  this  would  change  the  flow  field  that  should  be  analyzed  by 
the  boundary  layer  program.  Note:  Results  from  Case  4 in  Section  6.3 
suggest  that  the  influence  of  the  free  vortex  should  also  he  eliminated 
because  leading  edge  pressures  are  not  correct  until  vortex  converges  to 
correct  position.  Once  the  separation  point  is  established  though,  the 
vortex  sheet  and  the  secondary  vortex  would  improve  the  potential  flow 
pressure  distribution. 
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It  is  also  important  to  note  that  separation  did  not  occur  when 
turbulence  was  initiated  at  station  7 for  the  7.08  xIO5  Reynolds  number. 
(Note:  Transition  is  artificial  and  this  would  probably  not  cause 
turbulence.)  The  flow  was  attached  all  the  way  to  X/C  *.15  even  using  the 
attached  flow  pressure  distribution.  The  significance  of  this  is  that  at 
high  Reynolds  numbers  where  transition  occurs  naturally  there  will  be  no 
leading  edge  separation  and  therefore  no  leading  edge  vortex.  This  may  be 
a real  problem  for  configurations  like  the  space  shuttle  and  future  super- 
sonic transports  that  are  designed  to  take  advantage  of  the  leading  edge 
vortex.  Even  on  fighter  aircraft  where  the  wings  are  thin,  this  phenomenon 
will  determine  at  what  angle  the  vortex  forms  and  at  what  position  (affect- 
ing the  moments  and  stability).  It  is  very  important  that  a high  Reynolds 
number  boundary  layer  program  be  used  to  determine  where , how , and  how 
much  vorticity  (circulation)  is  fed  to  the  leading  edge  vortex.  (See 
Case  4,  Section  6.3.) 

h.h  Application  of  Boundary  Layer  Program 

As  noted  previously,  the  output  required  from  the  boundary  layer  is 
the  vorticity,  its  strength,  and  direction.  More  specifically,  the  rate 
of  circulation  shed  from  the  boundary  layer  in  all  three  directions  is  the 
required  output  from  the  boundary  layer  program. 

The  total  rate  of  circulation  for  each  component  of  vorticity  using 
boundary  layer  coordinates  transferred  to  box  coordinates  is  given  below. 
(Figure  4.1  illustrates  how  the  boundary  layer  thickness  occurs  along  the 
curved  coordinate.)  For  most  practical  purposes,  however,  6 can  be 
assumed  normal  to  the  surface  making 

dz  = dN  sin(9  - c)  and  dy  * dN  cos(9  - e) , 
where  9 is  the  airfoil  surface  slope.) 
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(The  subscripts  on  the  brackets  indicate  the  plane  to  which  the  vorticity 
is  normal.)  For  an  infinite  yawed  wing  dx  is  zero  over  the  boundary  layer 
height  leaving  only  one  term  in  each  of  the  last  two  brackets. 


If  the  direction  of  the  velocity  is  not  of  concern 
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The  circulation  per  unit  length,  r/l,  can  be  written  as 
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However,  the  product  of  ?/ L and  the  average  velocity  in  the  boundary  layer 
must  remain  constant  for  each  component  when  transferring  the  circulation 
to  the  viscous  vorticity  box.  That  is,  where  the  exact  position  of  the 
boundary  layer  cannot  be  maintained  because  of  the  large  grid  size  in  the 
box,  the  value  of  dF/dt  from  equation  (4.7)  must  be  maintained.  Therefore, 
the  velocity  at  the  nearest  grid  point  will  have  to  be  changed  to  the 
average  boundary  layer  velocity  and  equation  (4.8)  used  for  F/l.  Vorticity 
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per  unit  length  is  obtained  by  dividing  the  resulting  r/L  by  the  mesh 
length,  Ay,  of  the  box.  This  technique  eliminates  the  confusion  in 
transferring  vorticity  from  the  boundary  layer  to  the  box.  Section  5 
discusses  this  technique  further,  showing  the  relation  between  boundary 
layer  dimensions  and  box  dimensions. 

Originally  it  was  hoped  that  the  grid  points  of  the  box  would  actually 
fall  in  the  boundary  layer  as  shown  in  figure  4.1.  This  was  not  practical 
for  the  current  work  for  the  number  of  mesh  points  in  the  box  would  be  too 
large  for  practical  computer  time.  A grid  transformation  technique  would 
allow  closer  grid  points  near  the  boundary  layer  while  stretching  their 
distance  apart  further  away  from  the  boundary  layer.  Even  when  this  occurs 
though  it  will  be  necessary  to  apply  the  constant  rate  of  circulation 
concept  at  each  grid  point  which  must  contain  all  the  circulation  over  a 
finite  distance  of  the  boundary  layer. 
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5.  VISCOUS  BOX  CALCULATIONS 


5.1  Governing  Equations 

With  the  boundary  conditions  provided  by  the  potential  flow  and 
boundary  layer  calculations,  detailed  flow  features  are  computed  over 
finite  integration  domains  by  solving  the  Navi er-Stokes  equations.  The 
governing  equation  is  formulated  as  the  vorticity  transport  equation, 

^ • (Vu)  “ u ) • (vV)  + vV2w  (5-1) 

where  velocity  and  the  derived  variable  vorticity  given  by 

2 = V x V (5.2) 

are  used  as  the  primary  dependent  variables.  The  velocities  are  computed 
by  numerically  solving  the  equation 

V2V  - - V XU.  (5.3) 


If  the  physical  flow  phenomena  permits  the  assumption  that  diffusion 
is  negligible  in  one  of  the  coordinate  directions  (say  x-direction) 
compared  to  the  other  two  directions,  then  one  can  approximate  the  diffusion 
term  as 


This  renders  the  partial  differential  equations,  (5.1).  to  be  parabolic 
thus  enabling  the  usage  of  forward  marching  techniques.  The  integration 
domain  for  these  governing  equations  are  the  so-called  "Viscous  Boxes", 
the  orientations  of  which  are  shown  in  Figure  2.1.  The  details  of  the 
numerical  procedure  for  solving  these  governing  equations  may  be  found  in 
the  earlier  reports  (refs.  9,  10,  and  11).  In  the  following  sections 
refinements  and  further  modifications  introduced  in  the  solution  procedure 
are  discussed. 
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5.2  Initial  and  Boundary  Conditions 

The  parabolic  marching  technique  used  for  the  solution  of  the  vorticity 
transport  equation  requires  the  knowledge  of  the  initial  conditions  on 
vorticity  and  velocity  at  the  starting  face  of  Box  1 (see  fig.  2.1)  and  the 
boundary  conditions  on  the  rectangular  boundary  of  each  subsequent  cross- 
flow  plane  (see  fig.  2.2). 

The  following  procedure  is  employed  to  establish  the  initial  conditions. 
The  viscous  box  solution  is  started  at  the  first  face  of  the  first  box  with 
the  potential  flow  velocity  field  resolved  into  the  box  coordinates,  the 
boundary  layer  vorticity  that  feeds  into  the  box,  and  a lamb  vortex  of  low 
strength  placed  at  a location  compatible  with  the  potential  flow  solution. 
Holding  the  boundary  layer  input  of  vorticity  and  the  velocity  values  at  the 
boundary  constant,  the  solution  is  advanced  using  small  Ax  values  for  a few 
stations.  This  is  continued  until  (a)  a feeding  sheet  of  vorticity  connects 
the  boundary  layer  vorticity  to  the  Lamb  vortex,  and  (b)  the  total  vortex 
strength  becomes  equal  to  the  vortex  strength  employed  in  the  potential  flow 
model.  When  these  conditions  are  met,  the  solution  obtained  for  velocity 
and  vorticity  are  assigned  back  to  the  starting  face  as  the  required 
starting  conditions. 

As  was  pointed  out  in  Section  2,  the  extension  of  the  box  boundaries 
beyond  the  leading  edge  offers  possibilities  of  significantly  improving  the 
way  in  which  the  local  boundary  conditions  are  imposed  (see  figure  2.2). 

The  extended  region  of  the  no-slip  surface  permits  an  accurate  positioning 
of  the  boundary  layer  vorticity  input  points,  limited  only  by  the  grid  size 
of  box  program  in  relation  to  the  boundary  layer  thickness.  If  a cluster- 
ing transformation  technique  is  used,  the  cross-plane  geometry  would  even 
facilitate  the  boundary  layer  vorticity  profile  to  be  exact'y  incorporated 
at  the  box  boundary.  At  this  stage,  however,  no  grid  transformation  tech- 
nique is  used  and  the  vast  difference  between  the  grid  sizes  of  the  viscous 
box  and  the  boundary  layer  program  necessitates  certain  approximations. 
Essentially,  all  the  boundary  layer  vorticity  is  now  placed  at  one  point 
at  the  box  boundary  while  keeping  total  circulation  to  be  same  between  two 
cases.  The  details  of  interpreting  the  boundary  layer  results  compatible 
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with  the  box  grid  size  were  discussed  earlier  in  Section  4.  On  the  no-slip 
boundary  of  the  box,  velocity  is  specified  as  zero  and  the  vorticity  is 
determined  using  finite  differences.  The  vorticity  values  are  set  to  zero 
on  the  free  edges  of  the  box  and  the  velocities  are  obtained  by  resolving 
the  potential  flow  velocities  into  the  box  coordinates. 

5.3  Modifications  to  Finite  Difference  Technique 

Previously,  the  vorticity  box  program  employed  up-wind  differencing  for 
the  in-plane  derivatives  of  vorticity  while  solving  the  vorticity  transport 
equation,  (5.1).  The  advantages  (not  being  stability  limited  by  cell 
Reynolds  number,  being  able  to  preserve  "transport i ve  property")  and  dis- 
advantages (lower-order  formal  accuracy,  numerical  viscosity)  of  this 
difference  scheme  are  well  documented  in  open  literature  (see  for  example, 
ref.  18).  Whether  the  numerical  diffusion  introduced  by  the  upwind  dif- 
ferencing can  be  tolerated  or  not  depends  to  a large  extent  on  the  physics 
of  the  problem.  In  the  studies  of  vortex  growth  or  decay  in  a viscous 
medium,  errors  due  to  numerical  viscosity  become  critically  important. 

In  view  of  the  importance  of  understanding  the  errors  due  to  numerical 
viscosity,  an  extensive  study  of  error  analysis  was  undertaken  at  private 
expense  and  the  results  may  be  found  documented  in  reference  19.  This 
analysis  involves  the  study  of  a Lamb  vortex  positioned  axially  in  a viscous 
stream,  for  which  case  an  exact  analytical  solution  is  available.  As  shown 
in  figure  5.1.  the  fully  upwind  differencing  scheme  causes  the  vortex  to 
decay  very  rapidly  in  comparison  to  the  exact  solution.  It  was  also  found 
(ref.  19)  that  a differencing  scheme  like  central  differencing,  which  is  a 
higher-order  accurate  scheme  and  which  is  also  free  from  numerical  viscosity 
for  the  present  form  of  governing  equation,  gives  rise  to  "wiggles"  in  the 
solution  near  the  boundaries.  Thus,  a tradeoff  is  apparent  and  a combined 
central /upwind  differencing  is  found  to  give  the  best  result  (fig.  5.1). 
Under  the  new  scheme,  at  a given  grid  point  either  the  central  or  upwind 
differencing  is  used  depending  upon  the  value  of  a parameter,  $,  defined  by 
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where  C is  the  local  vorticity  value  and  Re  is  the  flow  Reynolds  number. 
Upon  studying  the  effects  of  various  grid  spacings  and  the  flow  Reynolds 
number  on  the  errors  for  the  ideal  Lamb  vortex  case,  the  optimum  value  of  $ 
has  been  found  to  be  in  the  range  25  s <f>  s 50.  If  one  chooses  a value  of 
<(>=30,  say,  then  the  numerical  scheme  is  to  employ  central  differencing  at 
grid  points  where  <(>>30  and  to  employ  upwind  differencing  at  points  where 
<(><30.  For  the  present  study  where  the  Reynolds  number  is  constant  and  the 
grid  spacing  at  box  face  is  also  a constant,  it  amounts  to  using  central 
differencing  where  the  vorticity  values  are  high  (like,  for  example,  near 
vortex  core  and  surface),  and  to  using  upwind  differencing  where  the 
vorticity  values  are  low  (typically,  the  region  near  free  boundaries). 

5. A Reconstruction  of  Free  Vortex  Geometry 

A further  refinement  introduced  in  the  viscous  box  program  compared  to 
the  prior  approach  is  in  the  method  of  analyzing  the  final  results  with  a 
view  to  reconstruct  the  free  vortex  geometry.  Previously,  all  the  vorticity 
in  each  cross-sectional  plane  was  integrated  to  obtain  a center  of  gravity 
position  for  the  vortex  and  to  determine  the  circulation  strength.  With 
such  an  approach,  it  is  not  possible  to  isolate  secondary  vortices,  if  any, 
from  the  primary  vortex.  (See  figure  6>.20  for  vorticity  contours.) 

In  the  new  version  of  the  box  program  contouring  techniques  have  been 
added  to  isolate  different  "vortex  regions."  Typically,  at  each  solution 
plane,  all  vorticity  contours  with  values  0+  and  0“  are  first  obtained, 
and  then  the  center  of  gravity  and  total  circulation  inside  each  such 
contour  are  determined.  After  repeating  this  procedure  for  each  cross- 
sectional  plane,  vortex  filament  geometry  and  variation  of  the  circulation 
strength  along  its  length  are  readily  reconstructed.  At  the  end  of  this 
procedure  the  vortex  geometry  definition  will  be  in  box  coordinates.  These 
are  then  transformed  to  Hess  coordinates  using  appropriate  transformation 
relation. 
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6.0  PROGRAM  INTERFACING  AND  SAMPLE  CALCULATIONS 

While  developing  the  integrated  computational  program,  some  valuable 
operating  experience  was  gained.  In  the  material  that  follows,  the  efforts 
are  applied  to  a 65°  delta  wing  as  described  in  Section  6.1.  The  efforts 
are  presented  in  cases  having  one  or  more  cycles.  A case  is  defined  as  a 
consistent  run  through  all  the  necessary  cycles  without  changing  the  method. 
A cycle  is  defined  as  one  complete  run  through  all  the  programs  as  shown  on 
the  chart  in  Section  2.0  with  the  end  being  pressures  and  forces  from  the 
potential  flow  model.  Four  cases  have  been  run  but  only  the  last  has  been 
consistent  through  each  of  its  three  cycles.  In  each  of  the  earlier  cases, 
the  operating  experience  suggested  a major  improvement  to  the  technique, 
and  it  was  deemed  appropriate  to  start  anew  rather  than  to  continue  with 
additional  cycles.  The  following  table  describes  the  cases  studied. 


TABLE  1 

OPERATING  EXPERIENCE 


Case 

Description 

Cycles 

Resul ts 

1 

Initial  estimate  of  free 
vortex  strength  high  by 
factor  2.35 

2 

Unrealistic  pressures 
and  vortex  growth 

2 

Initial  estimate  of 

Case  1 reduced  by  i 

1 

Vortex  growth  similar 
to  Case  1 ; Boundary 
layer  input  too  high 

3 

Boundary  layer  vorticity 
reduced  to  that  between 
6-6*;  Initial  cond  i- 
tions  same  as  Case  2 

1 

Vortex  growth  and 
position  qualitatively 
good;  Box  starting 
conditions  in  error 

1* 

Boundary  layer  rate  of 
circulation  input  re- 
duced by  rate  of  leaving 
the  box;  Starting 
conditions  of  box 
corrected;  Case  2 
initial  conditions 

3 

Boundary  layer  vorticity 
correct  only  for  potential 
flow  with  no  free  vortex; 

Must  converge  on  free 
vortex  position  and 
strength  before  changing 
boundary  layer;  Vortex 
position  must  be  repre- 
sented more  accurately 
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The  first  three  cases  are  discussed  in  Section  6.2.  They  are  presented 
in  this  report  primarily  to  illustrate  the  importance  of  some  of  the  parame- 
ters. Also  the  experience  gained  in  evolving  correct  operating  procedures 
will  assist  future  operators  of  the  program. 

Case  4,  which  incorporates  all  improvements  as  learned  from  the 
previous  case  is  discussed  in  Section  6.3. 

During  the  analysis  of  Cases  1 through  3,  it  was  found  that  an 
evaluation  of  the  cross-flow  velocity  (vx3)  at  the  free  vortex  axis  with 
the  potential  flow  program  would  enhance  convergence.  This  is  discussed 
i n Sect  ion  6.2.1. 

6.1  Model  Description 

In  order  to  develop  the  techniques  of  this  hybrid  potential/viscous 
flow  model,  a 65°  delta  wing  with  rounded  leading  edges  was  selected  for 
application.  This  wing  is  similar  to  one  tested  by  Henderson  (ref.  16) 
except  for  a small  trailing  edge  extension  on  his  wing.  The  reasons  for 
selecting  this  wing  are:  (1)  both  sharp  and  rounded  leading  edge  data  are 
available  for  overall  forces  and  (2)  the  wing  is  thin  enough  to  insure  that 
a leading  edge  vortex  is  formed  at  a relatively  low  angle  of  attack.  An 
angle  of  14°  was  selected  for  analysis  to  insure  a high  enough  angle  to  have 
vortex  formation  over  the  entire  leading  edge  but  not  too  high  to  cause 


bursting  over  the  wing. 

The  airfoil  selected  for  the  65°  delta  was  chosen  to  be  similar  to  that 
of  Henderson's  primarily  near  the  leading  edge.  The  streamwise  section  is  a 
NACA  0006  airfoil  reduced  to  a t/c  of  0.02302.  It  is  further  modified  from 
5%  of  x/c  to  the  leading  edge  to  correspond  to  the  largest  leading  edge 
radius  configuration  of  Henderson's  wing.  The  first  0 . 8 1 6%  of  the  leading 
edge  is  shown  in  figure  6.1.  From  0 . 8 1 6%  to  5.0$,  there  is  a straight  line 
fairing  into  the  reduced  thickness  0006  airfoil.  The  airfoil  is  maintained 
constant  from  root  to  tip,  and  since  the  tip  has  zero  dimensions  the 
outboard  stations  become  very  thin. 
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Figure  6.1.  Streamwise  Section  Leading  Edge  Shape  - t/c  = 0.02302 


These  thin  airfoils  made  it  necessary  to  have  the  viscous  box  configu- 
ration extended  beyond  the  leading  edge  as  shown  in  figures  6.2  through  6.5* 
The  face  and  end  of  each  box  as  shown  on  the  wing  in  figure  2.1  are 
illustrated.  The  box  layout  of  figure  2.1  is  used  for  the  analysis  and  the 
results  are  presented  in  this  section  of  the  report.  Box  1 is  not  used  in 
the  analysis  because  the  first  calculations  indicate  the  vortex  has  not 
formed,  i.e.  no  flow  separation. 

The  boxes  of  figures  6.2  through  6.5  are  shown  inclined  to  the  wing 
reference  plane  and  this  inclination,  e,  varies  with  each  position  of  the 
boxes.  It  was  necessary  to  do  this  in  order  to  have  the  bottom  of  the  box 
fit  the  airfoil  in  the  approximate  manner  shown.  Obviously,  a coordinate 
system  for  the  boxes  that  transforms  to  fit  the  airfoil  is  suggested.  The 
compromise  shown  here,  however,  should  not  significantly  affect  the  data 
since  the  airfoil  does  not  deviate  much  from  the  flat  box  floor.  In 
computing  the  boundary  velocities,  the  actual  box  positions  are  used; 
however,  the  viscous  box  calculations  ignore  the  small  twist  in  the  boxes 
caused  by  the  varying  e. 


6.2  Results  of  First  Three  Cases  of  Investigations 

In  Case  1 two  consistent  cycles  were  computed.  Initial  and  resulting 
vortex  strengths  and  positions  are  shown  in  figures  6.6  through  6.8.  The 
original  estimated  free  vortex  strength  and  position  is  shown  as  the  dashed 
lines  in  figures  6.6  and  6.8.  These  were  estimated  from  the  method  reported 
in  reference  11.  It  was  determined  only  after  two  cycles  of  the  program 
were  completed  that  the  original  estimate  of  vortex  strength  was  about  2.35 
times  its  correct  value.  Reasons  for  this  error  are  as  follows:  First, 
the  theoretically  estimated  values  for  vortex  strength  appear  to  be  the  sum 
for  both  sides  of  the  wing,  therefore,  a factor  of  two  is  determined.  Then, 
the  0.35  factor  is  an  adjustment  for  the  deviation  of  the  theory  from 
experiment.  Note:  This  estimating  technique  is  valid  only  for  a class  of 
sharp  edge  delta  wings,  and  further  corrections  are  needed  to  estimate  the 
effect  of  round  leading  edges  and  different  planfoms.  For  future  applica- 
tions a new  method,  using  the  circulation  from  the  boundary  layer,  is  being 
developed  to  obtain  the  initial  estimate  of  free  vortex  strength. 


As  noted  in  figures  6.6  and  6.8,  the  vortex  strength  and  vertical 
height  as  determined  from  the  first  and  second  cycles  differs  radically 
from  the  original  estimate.  This  behavior  was  originally  attributed  to  the 
highly  inaccurate  original  estimate  of  vortex  strength.  One  particularly 
significant  effect  was  the  pressure  distributions  used  to  predict  the 
leading  edge  boundary  layer  and  its  resulting  circulation.  Figure  6.9 
shows  the  erratic  circulation  from  the  boundary  layer  as  a function  of  the 
position  along  the  leading  edge.  Part  of  this  erratic  behavior  is  due  to 
the  pressures.  It  was  difficult  to  determine  from  the  pressure  distribution 
a proper  starting  position  (attachment  line  on  the  lower  surface).  The 
criteria  used  was  to  assume  that  the  point  where  the  most  positive  pressure 
occurred  was  on  the  attachment  line.  With  the  very  high  free  vortex 
strength,  the  lower  surface  pressure  coefficient  was  often  negative  making 
it  difficult  to  decide  where  the  attachment  occurs.  For  this  reason  and 
the  resulting  unusual  pressure  gradients,  the  boundary  layer  circulation 
prediction  was  erratic  as  a function  of  the  position  along  the  leading  edge. 
This  strange  behavior  may  also  be  due  in  part  to  the  leading  edge  changing 
from  very  rounded  near  the  apex  to  almost  sharp  near  the  wing  tip.  Separa- 
tion occurred  on  the  upper  surface  for  the  forward  part  of  the  wing  but  from 
about  50%  of  the  span  to  the  tip  the  flow  separates  at  the  leading  edge. 
Furthermore,  the  data  presented  in  figure  6.9  was  obtained  prior  to  the 
analysis  described  in  Section  A.O  for  the  22°  swept  wing.  In  that  case, 
more  attention  was  given  to  definition  of  points  near  the  separation  point. 

Figure  6.9  also  gives  the  boundary  layer  results  for  Case  2.  Here 
careful  attention  was  given  to  the  boundary  layer  analysis  and  the  pressures 
were  more  realistic.  Reasons  for  this  shape  are  given  subsequently  when 
comparing  with  Case  3. 

Instead  of  continuing  to  the  third  cycle,  a new  case  was  initiated 
with  a more  accurate  estimate  of  the  free  vortex  strength.  The  results 
of  the  first  cycle  of  this  Case  2 were  both  encouraging  and  disappointing. 
The  encouraging  point  is  that,  even  though  the  initial  vortex  strength 
estimate  was  reduced  by  a factor  of  1/2,  the  vortex  strength  and  position 
after  going  through  the  viscous  box  calculations  were  about  the  same  as  the 
first  cycle  of  Case  1.  This  would  say  that  the  viscous  box  results  are  not 


are  not  strongly  affected  by  the  initial  vortex  strength.  Of  course,  the 
discouraging  fact  was  that  the  box  results  were  just  as  radically  different 
from  the  original  estimated  values  as  that  for  Case  1.  It  was  then  necessary 
to  examine  other  problem  areas. 

After  some  review,  it  was  determined  that  an  improper  boundary  condition 
was  set  for  the  viscous  boxes.  Up  through  Case  3 the  shed  boundary  layer 
circulation,  T/L,  was  obtained  by  simply  averaging  the  vorticity  over  the 
boundary  layer  thickness  at  the  separation  point.  As  discussed  in  Section 
A. A,  the  correct  value  of  T/L  is  the  rate  of  circulation  divided  by  the 
average  boundary  layer  velocity.  Also,  up  through  Case  3,  no  vorticity  is 
allowed  to  leave  the  boxes.  This  condition  allows  too  much  of  the 
vorticity,  being  fed  into  the  box  from  the  boundary  layer,  to  stay  in  the 
box;  causing  the  free  vortex  to  be  too  strong.  The  most  accurate  means  to 
correct  this  problem  is  to  allow  vorticity  to  leave  the  box  at  its  trailing 
edge  near  the  airfoil  surface  as  boundary  layer  (see  figure  2.2).  One  means 
of  accomplishing  this  is  to  set  the  gradient  of  vorticity  along  the 
2-direction  of  the  box  to  zero;  i.e.  3ui/9z=0.  This  is  opposed  to  the 
current  practice  of  setting  5 = 0 at  the  edge  of  the  box.  This  concept  is 
equivalent  to  acknoqledging  that  part  of  the  vorticity  in  the  boundary  layer 
goes  into  forming  the  free  vortex  and  its  lift,  and  the  remaining  part  is 
that  which  corresponds  to  the  remaining  part  of  the  lift  created  by  the  wing 
surface.  In  viscous  flow  analysis,  the  wing  surface  vorticity  leaves  along 
the  entire  trailing  edge  to  join  with  the  free  vortex  vorticity  further 
downstream  of  the  wing.  This  fact  requires  that  some  vorticity  leave  the 
box.  If  one  knew  the  vorticity  leaving  the  box,  the  vorticity  from  the 
leading  edge  could  then  be  reduced  by  that  amount.  Then  the  condition  that 
tu  * 0 at  the  trailing  edge  of  the  box  would  be  a good  approximation.  Some 
attempts  at  reducing  the  leading  edge  vorticity  are  discussed. 

Case  3 was  initiated  with  reduced  vorticity  from  the  leading  edge  to 
simulate  the  effect  of  vorticity  leaving  the  box.  Here,  it  is  assumed  that 
the  portion  of  vorticity  going  into  the  box  should  be  that  above  the  dis- 
placement thickness,  i.e.  between  6 and  6*.  Using  the  initial  estimated 
free  vortex  strength  of  Case  2 (i  Case  1)  and  the  5 -6*  vorticity,  one  cycle 


of  the  program  was  run  as  Case  3.  The  ci rculation  or  vorticity  obtained 
from  the  boundary  layer  for  the  two  cases  is  presented  in  figure  6.10. 

Behavior  of  the  shedding  circulation  as  a function  of  X/CR  is  strange, 
but  informative.  Upon  investigation  it  was  determined  that  up  to  the  high 
point  at  box  4 face,  turbulence  occurs  before  separation.  The  remaining 
points  have  laminar  separation.  This  causes  some  of  the  strange  behavior, 
but  most  of  the  variation  in  circulation  as  shown  from  the  line  is  due 
to  the  relation  of  local  angle  of  attack  and  leading  edge  radius.  Circula- 
tion increases  along  the  leading  edge  initially  due  to  increasing  local 
angle  of  attack  and  reduction  in  the  leading  edge  radius.  Continuing  along 
the  leading  edge,  the  separation  moves  forward  to  an  optimum  position  from 
the  leading  edge,  (i.e.  peak  vorticity).  After  this,  the  radius  gets  so 
small  that  the  separation  rapidly  approaches  the  leading  edge  where  the 
vorticity  gets  smaller  because  turbulent  flow  has  not  developed  and  the 
limiting  condition  of  a sharp  leading  edge  has  been  obtained. 

The  information  gained  from  this  test  of  reduced  leading  edge  vorticity 
is  important.  The  level  and  shape  of  circulation  obtained  from  the  boundary 
layer  is  quite  different  for  Cases  2 and  3-  No  other  difference  occurs 
between  these  two  cases,  but  the  resulting  free  vortex  strength  and  position 
is  significantly  different  as  shown  in  figures  6.11  through  6.13.  In  each 
figure  the  results  look  qualitatively  better  for  Case  3*  The  only  real 
conclusion  though  is  that  the  boundary  layer  vorticity  is  an  extremely 
critical  parameter.  Careful  attention  must  be  given  to  computing  this 
parameter.  First,  the  airfoil  shape  and  pressure  distributions  must  be 
obtained  accurately,  and  second,  the  circulation  rate  should  be  used  to 
obtain  the  shed  leading  edge  circulation  (vorticity),  and  finally,  the 
boundary  conditions  on  the  viscous  box  must  provide  for  the  proper 
vorticity  leaving  the  box  as  well  as  entering  it. 

An  alternate  method  of  obtaining  the  rate  of  circulation  leaving  the 
box  is  used  for  Case  4 and  suggested  for  later  application.  It  is  assumed 
that  the  wing  surface  velocity  at  the  inboard  box  boundary,  as  calculated 
from  the  potential  flow  programs,  can  be  used  to  compute  the  existing  rate 
of  circulation  vector,  i.e. 
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Figure  6.11,  Free  Vortex  Vertical  Position  Cases  2 & 3 


Figure  6.13. 
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This  rate  could  be  used  to  compute  the  boundary  value  vorticity,  or  the  net 
rate  of  circulation  allowed  to  enter  the  viscous  box  could  then  be  the  dif- 
ference between  the  computed  boundary  layer  circulation  rate  and  the 
assumed  rate  of  equation  (6.1).  This,  of  course,  is  not  as  rigorous  as  the 
previously  suggested  method  which  sets  3ui/3z=0  at  the  box  inboard  side.  The 
latter  will  take  some  programming  changes  that  are  time  consuming,  but  the 
information  for  equation  (6.1)  is  already  available.  Case  4 is  computed 
using  equation  (6.1)  to  obtain  the  boundary  value  vorticity  leaving  the  box. 

6.2.1  Initial  Vortex  Position  and  Cross-Flow  Velocity 

In  the  computations  of  Cases  1 through  3,  a problem  associated  with 
vortex  position  was  discovered  only  after  going  through  Case  3-  This  was  a 
problem  causing  the  starting  conditions  at  the  face  of  box  2 to  be  incorrect 
at  the  face  of  box  2 to  be  incorrect.  At  this  face  there  is  provision  to 
allow  the  initial  free  vortex  to  be  input  in  its  three  orthogonal  directions. 
In  the  process  the  free  vortex  inclination  relative  to  the  wing  surface  was 
too  much.  This  caused  a large  distortion  of  the  flow  field  velocities  and 
errors  for  all  the  other  boxes.  This  problem  has  also  been  corrected  for 
Case  4. 

The  initial  estimate  of  vortex  position  can  be  improved  by  another 
technique,  which  evolved  during  the  program  interfacing  of  Cases  1 through 
3.  This  is  related  to  an  evaluation  of  the  cross-flow  velocity  at  the 
vortex  axis.  The  ultimate  result  of  convergence  with  this  program  is  no 
force  on  the  vortex  and  no  change  in  vortex  strength  between  cycles.  In 
mathematical  terms  no  force  on  the  vortex  requires: 

V x io  3 0,  at  vortex  axis.  (6.2) 

Equation  (6.2)  is  not  satisfied  until  final  convergence  of  the  cycles; 
therefore,  neither  the  pressures  nor  forces  calculated  by  the  potential  flow 
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model  are  correct  until  convergence.  This  is  shown  by  examination  of  the 
integral  equation  of  motion  for  steady  flow  and  negligible  body  forces; 


V 


x 


(6.3) 


With  circulation  (free  vortex)  in  the  flow  field  having  a net  force 
(VxW^O),  the  Bernoulli  equation  for  pressures  in  irrotational  flow  is  no 
longer  valid.  This  causes  the  preconvergent  pressure  integration  technique 
of  the  AIR  program  to  give  erroneous  forces. 

As  was  noted  previously,  if  very  large  errors  are  present  in  the 
estimated  free  vortex  strength  and  position,  the  pressures  and  interference 
effects  between  wing  and  vortex  are  largely  in  error.  To  avoid  this  problem 
before  going  through  all  the  viscous  analyses,  a check  on  V xiL  should  be 
made  with  the  potential  flow  program  and  the  vortex  position  changed 
accordingly.  Figure  6.14  illustrates  such  an  analysis.  The  velocity 
normal  to  the  vortex  axis  and  parallel  to  the  wing  reference  plane  is  given 
as  a function  of  the  X/Cr.  Results  from  the  initial  position  of  vortex  in 
Case  1 indicate  the  vortex  wants  to  move  vertically  from  the  surface  near 
the  apex  and  towards  the  surface  as  the  vortex  moves  to  the  wing  trailing 
edge.  If  this  tendency  had  been  relieved  by  moving  the  vortex  in  its 
desired  direction  prior  to  going  through  the  viscous  analyses,  the  new 
starting  conditions  for  the  viscous  analyses  would  have  been  more  compatible. 
This  was  not  done  and  figure  6.6  shows  the  viscous  box  results  did  move  the 
vortex  in  the  directions  indicated  by  the  Vxu  analysis.  However,  because 
of  the  extremely  oversized  estimate  of  vortex  strength  and  the  other 
inaccurate  conditions  mentioned  in  Section  6.2,  the  viscous  box  overshot 
the  positions  required  by  the  potential  flow  model.  The  results  of  the 
potential  flow  V xu  analysis  for  the  end  of  cycle  1 are  also  shown  in 
figure  6.14.  The  increase  in  vortex  height  near  the  apex  did  relieve  the 
force  on  the  vortex  some  but  the  high  vortex  strength  produced  by  the  box 
results  at  the  outboard  part  of  the  wing  caused  even  more  discrepancy  in 
the  vortex  position.  A part  of  this  discrepancy  in  the  box  3 and  4 results 
is  due  to  inaccurate  starting  conditions  in  box  2.  The  cycle  2 vortex 
position  shown  in  figure  6.6  shows  an  encouraging  move  of  the  vortex  back 
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towards  the  initial  position  in  the  box  k area.  Unfortunately  with  so  many 
problems  this  case  was  instructive  but  would  probably  never  converge. 


The  check  on  V xu  is  applied  after  potential  flow  runs  in  Case  4. 

This  loop  in  the  program  hopefully  would  have  the  effect  of  approximating  a 
more  accurate  vortex  position  for  a constant  vortex  strength.  The  strength 
and  position  is  then  varied  by  the  viscous  box  portion  of  the  cycle. 
Unfortunately  in  Case  **  the  vortex  was  moved  in  the  wrong  direction.  Further 
investigation  is  required  before  this  technique  can  be  used  successfully. 

6.3  Case  4 - Incorporation  of  Program  Corrections 
Developed  in  Cases  1 through  3 

The  techniques  learned  in  Cases  1 through  3 and  applied  in  Case  **  are 
as  follows: 

o Apply  improved  original  estimate  of  free  vortex  strength. 

o Check  for  0 at  the  free  vortex  after  each  potential  flow 

run  to  improve  convergence. 

o Maintain  turbulent  boundary  layer  separation  along  entire  leading 
edge. 

o Maintain  the  rate  of  circulation  coming  from  the  leading  edge 
rather  than  the  absolute  value  of  circulation  per  unit  length. 

o Compute  the  rate  of  circulation  leaving  the  box  from  equation  (6.1) 
and  apply  as  boundary  condition  at  inboard  side  of  each  box. 

o Use  average  rate  of  circulation  between  two  box  X stations  to 
account  for  boundary  vorticity  (circulation)  entering  between 
box  stations. 

o Apply  correct  vortex  inclination  at  starting  position  of  box  2. 

Each  of  the  above  techniques  were  successfully  applied  to  Case  A. 

Also,  improvements  were  made  in  program  operating  procedures.  All  manual 
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operations  were  eliminated  except  evaluation  of  the  boundary  layer  output 
and  its  input  to  the  viscous  box  program.  All  potential  flow  runs  were  put 
into  one  run  stream;  that  is,  the  Hess  and  vortex  lattice  programs  were 
combined  so  that  with  an  input  of  free  vortex  strength  and  position  the  out- 
put gives  velocities  at  the  box  positions,  pressures  and  leading  edge  for 
boundary  layer,  total  force  and  moment  coefficients,  and  surface  pressures. 
This  output  is  filed  on  the  computer  for  access  by  the  boundary  layer  and 
viscous  box  program.  The  potential  flow  program  also  checks  the  velocities 
at  the  free  vortex  axis  and  makes  an  adjustment  in  position  of  the  vortex 
axis  according  to  the  cross-flow  velocities,  V xw.  This  adjustment  proved 
to  be  in  error  and  needs  further  investigation  before  efforts  subsequent 
to  Case  A. 

Results  from  Case  4 provided  considerable  insight  to  convergence 
techniques.  Although  Case  A is  not  convergent  in  three  cycles,  the  effort 
has  provided  the  knowledge  necessary  to  obtain  convergence.  The  most  signi- 
ficant error  made  in  Case  4 is  still  in  the  application  of  the  boundary 
layer  output.  The  fact  that  surface  pressures  are  not  correct  when  Vxdi^O 
cannot  be  ignored.  The  only  conditions  where  the  boundary  output  is 
correct  is  for  the  potential  flow  without  a free  vortex  and  for  Vxto  = 0 in 
the  final  converged  cycle.  As  can  be  seen  in  figure  6.15,  the  boundary 
layer  circulation  shed  into  the  viscous  box  increases  for  each  cycle.  This 
causes  the  total  free  vortex  strength,  figure  6.16  (as  predicted  by  the 
viscous  box),  to  increase  with  each  cycle;  i.e.  a divergent  condition.  The 
only  reasons  for  circulation  from  the  boundary  layer  to  increase  with  each 
cycle  are:  (1)  the  leading  edge  pressure  increased  due  to  Vxu^O  and  (2) 
probably  the  free  vortex  sheet,  and  secondary  vortex  need  better  definition 
in  the  potential  flow  model. 

It  is  quite  apparent  now  that  the  solution  to  the  above  problem  is 
simply  to  relax  the  boundary  layer  input;  that  is  maintain  the  boundary 
layer  input  from  the  potential  flow  without  free  vortex  for  all  cycles  until 
the  vortex  position  and  strength  converges.  Then  the  vortex  sheet  and 
secondary  vortex  can  be  defined  .more  accurately  resulting  in  accurate 
leading  edge  pressures  which  can  be  used  in  the  boundary  layer  program  to 
refine  the  shed  vorticity.  Current  experience  suggests  that  convergence 
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CYCLE  3 


Figure  6.16  Vortex  Strength  Output  From  Viscous  Box 


will  occur  very  rapidly  for  the  first  phase  of  above.  In  fact,  it  is 
conceivable  that  no  manual  interaction  will  be  required  in  the  computer  as 
long  as  the  boundary  layer  input  is  constant. 

In  spite  of  the  lack  of  convergence.  Case  4 did  indicate  very  realistic 
output.  When  one  considers  the  convection  direction  for  vorticity  shedding 
from  the  leading  edge,  figure  6.17,  and  the  vortex  position  in  figure  6.18, 
there  is  a realistic  tendency  for  the  vortex  y-position  to  turn  inward 
towards  the  root  and  z tends  to  go  to  zero  near  the  beginning  of  box  2 (see 
line  marked  improved  estimate).  At  the  face  of  box  2 the  convection  direc- 
tion (obtained  from  boundary  layer  output)  is  nearly  parallel  to  the  airfoil 
all  the  way  to  the  root.  Therefore,  very  little  vortex  height  should  be 
expected  and  then  it  should  be  near  the  root.  Unfortunately,  in  all  three 
cycles  of  Case  4,  the  vortex  was  restrained  to  the  initial  conditions  at 
the  face  of  box  2.  It  should  have  been  allowed  to  relax  to  the  "improved 
estimate"  even  for  cycle  2.  This  could  not  easily  be  accomplished  with  the 
geometric  restraints  in  the  program. 

The  current  program  did  not  allow  the  vortex  position  to  be  defined 
accurately  enough  . Even  though  the  vortex  strength  could  be  defined  in 
twenty  segments  the  vortex  position  could  be  defined  in  only  four  linear 
segments.  This  is  an  unnecessary  restraint  which  can  be  removed.  This 
restraint  was  also  a problem  when  testing  for  Vx£3=0.  It  was  found  here 
that  the  overall  lift  coefficient  and  surface  pressures  were  very  sensitive 
to  the  vortex  position  and  the  desired  accuracy  could  not  be  obtained  with 
only  four  segments. 

Another  very  realistic  result  is  indicated  with  Case  4.  The  convection 
direction,  9,  of  the  boundary  layer  vorticity  is  given  for  all  cycles  of 
Case  4 in  figure  6.17.  All  of  these  boundary  layer  outputs  are  computed  for 
a Rn«85.6x1  06.  Also  shown  is  the  convection  directions  for  R^  » 85.6x10® 
and  8.56  xIO6  for  the  boundary  layer  from  the  potential  flow  without  a free 
vortex.  All  of  these  show  large  differences  in  convection  direction.  When 
comparing  the  dashed  line  at  R^  ■ 85. 6 x 1 06  with  those  of  Case  4 at  the  same 
R^,  it  is  easy  to  see  that  the  large  values  of  9 come  from  increasing  the 
free  vortex  strength  and  this  is  not  accurate  as  discussed  earlier. 
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However,  a large  9 appears  necessary  for  any  significant  vortex  formation. 
The  realistic  conclusion  here  is  that  the  high  Rfj  of  85.6  xlO6  delays  large 
vortex  formation  to  higher  angle  of  attack.  On  the  other  hand  the  data  for 
Rn  = 8.56  xlO6  show  that  9 is  probably  high  enough  for  vortex  formation  at 
this  a = l4°.  A Case  5 is  now  being  computed  to  test  this  low  Rn  case 
including  all  the  other  corrections  mentioned  above.  A Case  6 is  also 
recommended  using  the  high  R^  boundary  layer  without  free  vortex  but  with 
smaller  viscous  boxes  or  finer  grids  (much  smaller  vortex  is  expected). 

An  example  of  the  effect  of  large  9 is  illustrated  for  cycle  3 in  the 
flow  field  plots  of  figure  6.19.  Only  in  box  3 does  the  flow  fie  d actually 
reverse  below  the  vortex.  This  is  because  9 is  large  or  in  other  *ords  the 
W component  of  boundary  layer  convection  velocity  is  small.  In  bo < 4 some 
of  the  detail  is  lost  due  to  large  grid  size  as  well  as  0 becoming  small 
again.  The  vorticity  contours  of  figure  6.20  show  some  secondary  vorticity 
which  should  be  represented  in  potential  flow. 

In  conclusion  the  analysis  of  Case  4 has  provided  the  following 
d i rect ion : 

(1)  Relax  boundary  layer  input  until  vortex  strength  and 
position  converge.  Then  check  boundary  layer  with 
accurate  vortex  system  represented. 

(2)  Improve  accuracy  of  vortex  position  and  represent 
vortex  feeding  sheet  and  secondary  vortex  in 
potential  flow. 

(3)  Carefully  consider  effect  of  Reynolds  number  in 
estimating  initial  vortex  position  and  strength. 
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7.0  APPLICATIONS 


7.1  General  Capability 

The  technique  developed  and  currently  discussed  in  this  report  is 
considered  to  be  an  optimum  combination  of  state-of-the-art  methods  which 
in  their  combination  allow  realistic  evaluation  and  design  of  the  critical 
geometry  for  vortex  control.  Vortex  control  can  then  be  used  to  improve 
high  angle  maneuvers  by  both  overall  lift  increase  and  aircraft  control. 

It  is  further  believed  that  through  optimized  design  that  a powerful  leading 
edge  vortex  can  be  obtained  at  the  cost  of  only  a small  fraction  of  the 
leading  edge  thrust.  The  hybrid  technique  discussed  here  is  developed  for 
arbitrary  leading  edge  shapes,  not  just  sharp  edges.  This  allows  for 
optimized  variation  of  leading  edge  shape  across  the  span  either  for 
roll/pitch  control  or  lift-to-drag  ratio.  It  is  not  difficult  to  envision 
a variable  geometry  leading  edge  (such  as  proposed  in  Supersonic  Commercial 
Aircraft  Research^®,  SCAR)  that  will  provide  optimization  for  both  vortex 
1 i ft  and  L/D. 

The  technique  developed  here  is  not  limited  by  scale  effects  even 
though  computational  differencing  techniques  are  used  in  the  viscous  vortex 
analysis.  All  of  the  boundary  conditions  to  the  viscous  box  are  available 
from  large-scale  techniques,  i.e.  both  the  potential  flow  and  boundary 
layer  programs.  In  this  sense  there  is  a practical  combination.  If  the 
viscous  box  was  modified  to  account  for  boundary  layer  around  the  leading 
edge,  the  Reynolds  number  would  be  limited  to  small  scale.  It  is  further 
believed  that  scale  effects  within  the  viscous  box  are  not  seriously  limit- 
ing the  application.  In  other  words  the  vortex  strength  and  position  is 
more  dependent  on  the  boundary  conditions  than  the  internal  flow  scale 
effects  or  grid  size.  In  contrast  though  the  viscous  effects  within  the 
box  are  necessary  to  define  both  the  primary  and  secondary  vortex  without 
resorting  to  various  assumptions  or  empiricism.  Currently  turbulence  is 
not  considered  in  the  viscous  box  analysis,  and  this  may  be  necessary  to 
consider  vortex  bursting. 

Vortex  bursting  cannot  currently  be  predicted  with  the  technique 
discussed  here,  but  there  is  good  prospect  that  the  parameters  are  available 
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from  the  method  to  use  in  empirical  vortex  bursting  analyses.  One  example 
is  the  swirl  parameter  discovered  by  Wilson^'.  If  incipient  bursting  does 
not  occur  in  the  viscous  box  before  the  critical  swirl  parameter  is  reached, 
then  the  true  vortex  burst  may  be  predicted  by  use  of  Wilson's  parameter. 
Incipient  burst  is  when  reverse  flow  occurs  in  the  direction  of  the  vortex 
axis.  It  is  here  that  turbulence  may  be  necessary  to  realistically  delay 
the  incipient  burst. 

Even  though  there  are  several  programs  involved  in  the  computational 
procedures  of  this  report,  the  computer  time  is  relatively  small.  Currently, 
it  costs  $82  or  less  than  5 minutes  of  NASA-Ames  CDC  7600  computer  time  to 
run  through  one  cycle  of  the  procedure;  i.e.  all  programs.  It  is  estimated 
that,  with  experience  in  operating  technique,  nominally  only  three  cycles 
will  be  required  to  converge  the  solution  with  six  cycles  being  the  most 
required.  This  is  15  to  20  minutes  of  computer  time  which  is  considerably 
cheaper  than  the  fabrication  and  modification  of  wind  tunnel  models  required 
to  get  only  a part  of  the  information  given  by  the  computer  method.  This  is 
not  to  say  the  computer  method  will  eliminate  wind  tunnel  test,  but  it  will 
allow  a more  optimized  wind  tunnel  model  to  be  built.  An  alternate  use  of 
the  method  is  to  extrapolate  results  of  the  wind  tunnel  tests  to  large 
scale.  At  the  same  time,  synthesis  of  the  various  parametric  effects  will 
provide  knowledge  of  the  details  giving  a better  understanding  of  the 
mechanisms  and  means  toward  design  finesse. 

Details  of  the  flow  field  will  be  seen  even  more  accurately  in  the 
future  with  further  development  of  the  method  of  this  report.  For  this 
method  is  in  line  to  take  advantage  of  the  future,  rapidly  growing,  compu- 
tational techniques.  As  will  be  discussed  in  the  recommendations,  there  are 
already  several  improvements  that  will  make  the  method  applicable  to  more 
configurations  and  make  it  more  accurate.  It  is  believed  that  even  now 
application  can  currently  be  made  to  many  configurational  problems,  and 
this  will  improve  with  update  of  the  method. 

Some  of  the  practical  applications  suggested  are  as  follows: 
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(1)  High  Reynolds  number  optimum  spanwise  variation  in  leading  edge 
design  to  maximize  vortex  lift  to  drag  ratio  for  both  fighters  and 
supersonic  transports. 

(2)  Optimized  variable  shape  leading  edge  for  high  angie-of-attack  roll 
and  pi tch  control . 

(3)  Canard  vortex  wing  interference  design. 

(k)  Analyses  of  leading  edge  vortex  effect  on  wings  at  high  angles  of 
attack  i n s ides  lip. 

(5)  Analyses  of  empennage  surfaces  for  vortex  flow  effects  while  in 
sidesl ip. 

(6)  Design  of  optimized  large  leading  edge  radius  for  configurations  such 
as  the  Space  Shuttle  for  optimum  vortex  lift. 

(7)  Analyses  and  optimization  of  vortex  formation  on  fighter  aircraft  with 
swept  forward  wings  and  subsequent  effects  on  their  empennage. 
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8.0  CONCLUSIONS 


o A theoretical  technique  has  been  developed  through  the  hybrid  combination 
of  potential  and  viscous  flow  models  to  predict  and  analyze  the  flow 
field  and  forces  due  to  a discrete  leading  edge  free  vortex  system.  This 
vortex  system  can  be  formed  from  arbitrary  wing  planforms  with  arbitrary 
edge  shapes. 

o Some  experience  has  been  gained  and  presented  herein  of  the  operating 
techniques.  These  techniques  are  still  evolving.  Experience  has  led  to 
improvements  of  the  techniques  and  to  recommendations  for  further 
i mprovement . 

o Current  applications  of  the  technique  have  been  applied  to  a 65°  delta 
wing  with  finite  thickness  and  spanwise  variation  in  leading  edge  radius. 
Absolute  convergence  of  the  technique  has  not  been  obtained,  but  the  . 
current  success  indicates  good  probability  of  convergence.  Additional 
cycles  and  experience  with  the  operation  will  provide  improved 
convergence  techniques. 

o Results  to  date  show  that  knowledge  of  rate  and  direction  of  leading  edge 
vorticity  is  an  absolute  necessity.  The  overall  vortex  lift  and  leading 
edge  thrust  are  very  sensitive  to  this  parameter,  and  this  parameter  is 
highly  sensitive  to  leading  edge  shape,  local  angle  of  attack,  Reynolds 
number  and  position  of  transition  from  laminar  to  turbulent  flow. 
Turbulent  flow  separation  provides  higher  vorticity  and  consequently 
higher  vortex  lift. 

o The  methods  used  herein  do  not  require  any  empiricism  except  probably 
for  future  application  to  vortex  burst  prediction. 

o Currently  the  methods  are  best  handled  with  engineering  interpretation 
between  some  of  the  steps;  however,  future  operating  experience  should 
eliminate  this  necessity. 
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9.0  RECOMMENDATIONS  FOR  FURTHER  EFFORT 


Recommendations  for  improvement  fall  into  two  categories.  One  makes 
significant  program  changes  to  current  programs  and  the  other  makes  the 
current  techniques  unified  as  well  as  conducting  applications  to  determine 
better  operating  techniques  and  make  sensitivity  studies. 

In  the  first  category,  six  (6)  program  changes  are  recommended. 

These  are: 

(1)  Change  coordinate  system  of  viscous  box  from  rectangular  to 
spherical  to  eliminate  jumping  in  box  size  and  to  compute  more  accurately 
in  a coordinate  system  natural  to  vortex  growth. 

(2)  Develop  a grid  transformation  technique  with  new  coordinates  to 
allow  more  detail  near  the  surface  of  the  wing,  improving  secondary  vortex 
definition  and  boundary  layer  input. 

(3)  Set  the  boundary  condition  on  the  inboard  grid  boundary  to  allow 
vorticity  to  exit;  i.e.  set  the  gradient  of  vorticity  to  zero  instead  of 
the  vorticity  to  zero  at  the  exit  boundary. 

(4)  Develop  a good  representation  of  vortex  feeding  sheet  for  the 
potential  flow  model  using  the  results  from  the  viscous  flow  analysis. 

This  will  require  a distinction  between  the  vorticity  in  the  discrete 
vortex  and  that  in  the  sheet.  Defining  the  velocities  in  cylindrical 
coordinates,  using  the  point  where  both  horizontal  and  vertical  components 
of  velocity  go  to  zero  as  the  coordinate  center,  will  help  define  the 
vortex  core  and  its  circulation.  This  will  also  provide  the  swirl 
parameter  needed  for  vortex  burst  analysis. 

(5)  In  the  potential  flow  vortex  lattice  model,  set  up  panels  in  the 
vortex  sheet  to  provide  no  flow-through  conditions. 

(6)  In  the  boundary  layer  program,  provide  means  for  automatically 
transitioning  from  laminar  to  turbulent  flow  using  the  best  physics  possible 
at  this  time. 
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In  the  second  category  a task  should  be  set  up  to  take  advantage  of 
the  current  programs.  The  following  is  recommended: 

(1)  Combine  all  programs  into  one,  requiring  as  little  manual  inter- 
action as  possible.  Where  manual  interaction  is  required,  use  computer 
graphics. 

(2)  Incorporate  grid  transformation  for  rectangular  coordinates  for 
same  purpose  as  (2)  above. 

(3)  Conduct  sensitivity  analysis  of  box  size,  Reynolds  number,  and 
lower  surface  attachment  line. 

(4)  Increase  number  of  potential  flow  model  panels  to  improve  pressure 
distribution  around  leading  edge. 

(5)  Fabricate  model  with  leading  edge  radius  such  as  the  one  in  this 
report  and  investigate  flow  field  with  laser  velocimeter  to  compare  results 
with  theory.  Test  at  an  angle  of  attack  without  vortex  burst  and  one  with 
burst.  Analyze  with  theory  for  both  cases. 


(6)  Investigate  methods  for  predicting  vortex  burst  using  theory 
and  empiricism  as  discussed  in  the  text. 
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